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PERIODIC  ORBITS  ABOUT  ROTATING  ASTEROIDS  IN  FREE  SPACE 


I .  Introduction 

Thousands  of  asteroids  revolve  around  the  Sun,  mainly 
between  the  orbits  of  Mars  and  Jupiter.  Some  of  these  are 
large  enough  to  consider  orbiting  with  a  small  probe  or 
satellite.  There  is  also  evidence  that  at  least  a  few 
asteroids  have  their  own  natural  satellites  (1:9-10).  For 
these  reasons,  this  study  investigates  families  of  orbits 
about  such  bodies. 

As  is  common  practice,  it  is  assumed  that  asteroids 
have  a  triaxial  ellipsoid  shape  with  axes  a  >  b  >  c  and  a 
rotation  about  the  shortest  axis  (7:436).  (See  Figure  1.) 
Further,  the  g r a v i t a t i on a  1  perturbations  of  the  Sun, 
planets,  and  other  asteroids  are  assumed  to  be  negligible 
(free  space).  Orbits  around  such  bodies  have  been  investi¬ 
gated  in  the  past,  but  no  previous  attempt  has  been  made  to 
find  the  equations  of  motion  in  such  a  way  as  to  be  condu¬ 
cive  to  numerical  solution  (3:37-38;  8:707-710;  9:75-84). 
Also,  because  these  earlier  studies  did  not  compute  actual 
orbits,  they  made  no  attempt  to  describe  the  appearance  of 
such  orbits. 


1 


Figure  1:  Asteroid  Geometry 


The  specific  task  in  this  study  is  to  derive  the  equa¬ 
tions  of  motion  in  such  a  way  as  to  facilitate  numerical 
solution  for  orbits  about  an  arbitrary  triaxial  ellipsoid  in 
free  space.  In  process  of  finding  these  orbits,  the  equi¬ 
librium  points  for  such  bodies  will  also  be  investigated. 
Finally,  as  an  example,  the  results  for  a  fictitious,  but 
realistic,  asteroid  will  be  presented  in  a  graphical  format. 
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II. 


Problem  Dynamics 
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The  problem  to  be  solved  is  that  of  finding  the  motion 
of  a  small  satellite  orbiting  an  asteroid.  It  is  assumed 
that  the  asteroid  can  be  modeled  as  a  triaxial  ellipsoid, 
rotating  about  its  shortest  axis  with  some  angular  velocity. 
If  the  satellite  mass  is  so  small  as  to  not  effect  the 
motion  of  the  asteroid,  then  it  can  be  assumed  that  the 
center  of  mass  for  the  system  is  fixed  in  inertial  space  and 
is  located  at  the  center  of  mass  of  the  asteroid.  Further, 
the  problem  is  restricted  to  asteroids  in  free  space. 

Equations  of  Motion 

The  equations  of  motion  are  to  be  numerically  solved, 
so  they  are  derived  with  this  in  mind.  Most  numerical 
integration  packages  are  written  to  solve  first  order 
differential  equations;  therefore,  the  Hamilton  canonical 


V, 


equations  are  derived  and  used  to  describe  the  motion. 

The  first  step  in  finding  the  equations  of  motion  is  to 
form  the  Hamiltonian 


H  =£Pi^i  "  L 


(2-1) 


where  p^  and  q^  are  the  generalized  momenta  and  velocities, 
respectively.  The  last  term,  L,  denotes  the  Lagrangian, 
defined  simply  as 


L  =  T  -  V 


(2-21 


Figure  2:  Gravitational  Potential 


with  T  and  V  being  the  kinetic  energy  and  potential  energy 
of  the  satellite.  It  is  these  two  terms  that  must  be 
determined  before  the  Hamiltonian  can  be  assembled. 

Potential  Energy.  The  gravitational  attraction  of  the 
asteroid  on  the  satellite  provides  the  potential  energy  and 
an  expression  for  this  can  be  found  (4:431,  434-436). 
Figure  2  shows  the  body-fixed  (rotating)  coordinate  system 
to  be  used  for  this  calculation.  From  the  figure,  it  is 
clear  this  can  be  written  in  the  integral  form 


V 


/' 


-Gm-,  /  (dm-,)  /r 


(2-3) 


where  r  =  |  K  -  J5 1  is  the  distance  shown  in  Figure  2  and  G 
is  the  Universal  Gravitational  constant. 


4 


The  1/r  term  in  Eq.  (2-3)  can  be  written: 


r'1  .  r’1 


2R  •  p 


(0 


-1/2 


(2-4) 


If  it  is  assumed  that  the  asteroid  dimensions  are  smaller 


than  the  distance  between  the  mass  centers,  then  the  right- 


hand  side  can  be  expanded  in  a  power  series: 


-1  -1 
r  =  R 


R  •  p  1  /p' 


2  \R 


3  2R  •  p  /p 


-  -  l-l  +  - 


8  R‘ 


•0 


+  •  •  • 


(2-5) 


Carrying  out  the  algebra,  this  can  be  rewritten  as: 


-1  -1 
r  =  R  + 


R.p  1  p2  3  (R  •  p)  2 


R3  2  R3 


+  order 


2  R- 


(0) 


(2-6) 


Using  the  body-fixed  reference  frame  previously  shown 


in  Figure  2  and  applying  the  following  definitions 


p  =  xex  +  yey  +  zez 


R  =  R  (le  +me  +ne) 
x  y  z 


(2-7) 


R  •  p  =  R  (lx  +  my  +  nz) 


where  1,  m,  and  n  are  the  appropriate  direction  cosines,  Eq. 


(2-6)  can  be  substituted  into  Eq.  (2-3)  to  yield  an 


>  *» 
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expression  for  the  potential.  After  ignoring  the  higher 
order  terms  and  some  simplification,  this  is  written  as: 


“Gml  f  Gml  f 

V  =  -  /  dm  -  — — —  /  (lx  +  my  +  nz)  dm, 

r  J  2  R2  J  <. 


Gml  f  2 

— -  /  [  3  (lx  +  my  +nz) 

2R3  J 


2  2  2 
(x  +  y  +  z  )  ]  dm. 


(2-8) 


The  first  integral  is  trivial: 

-Gmi  y.  -Gmim2 

-  /  drn  =  - 


(2-9) 


The  second  integral  becomes  zero  when  the  origin  of  the 
axis  system  is  placed  at  the  center  of  mass. 

The  final  term  of  Eq.  (2-8)  can  be  expanded  into 
easily  recognizable  quantities: 


-  f  [  3  (lx  +  my  +  ny)  2  -  (x2  +  y2  +  z2)  ]  dm_ 

2R3  J  ^ 


Gml  T  2  f  2  2  f  2 

— -  (31  -  1)  J  x  dm2  +  ( 3rn  -  1)J  yz  dm2 

n2  -  1) J'  z2  dm2  +  61m  J' xy  dm2 
n  J * xz  dm2  +  6mn  J' yz  dm2j 


+  (3n2  -  1)  j  z2  dm2  +  61m  j xy  dm 


+  61n  J  xz  dm2  +  6mn  J  yz  dm 


(2-10) 


This  can  be  simplified  by  noting 


V 
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dm  2 


dm2 


1/2  y  ( (x2  +  z2)  +  (x2  +  y2)  -  (y2  +  z2)]  dm2 

!/2  dyy  +  !zz  + 

1/2 /  [(x2  +  y2)  +  (y2  +  z2)  -  (x2  +  z2)]  dm2 
^  2  ^zz  +  Ixx  +  Iyy ^ 

1/2  J'  [ < y 2  +  z2)  +  (x2  +  z2)  -  (x2  +  y2)]  dm2 
l^2  ^xx  +  *yy  +  *zz^ 


/ 

/ 

/ 


xy  dm 


xz  dm 


yz  dm 


2 


2 


2 


xy 


xz 


yz 


(2-11) 


where  Ixx»  Iyy,  Izz  are  the  mass  moments  of  inertia  and  Ixy, 
Ixz,  IyZ  are  the  mass  products  of  inertia  (4:435-436).  If 
it  is  assumed  that,  in  addition  to  originating  at  the  center 
of  mass,  the  axes  are  arranged  such  that  Ixx,  Iyy,  and  Izz 
are  principal  moments  of  inertia,  then  Ixy  =  Ixz  =  Iyz  =  0. 


Clearly,  Eq.  (2-8)  can  now  be  written  as: 


-Gm^  rn  ^ 

Gm^ 

2 

=  - 

— 

-  [ 

(  3 1 z  - 

n 

R 

4R3 

yy 

+  (3m2  - 

1) 

^xx 

+  *zz 

■  Iyy) 

+  ( 3n2  - 

1) 

^  *xx 

+  xyy 

Izz^  ^ 

I 

zz 


I 

XX 


) 


(2-12a) 


The  body-fixed  spherical  coordinates  shown  in  Figure  3 
will  prove  to  be  more  convenient,  so  Eq.  (2-12a)  is 
rewritten 


-Guii  Bin  Gm,  _  - 

V  - - -  [  (3cos>  cos  \  -  1)  (I  +  I  -  I  ) 

R  4R3  yy  ZZ  xx 

+  ( 3cos20  sin2\  -  1)  (Ixx  +  Izz  -  I  ) 

+  ( 3s in20  -  1)  (Ixx  +  Iyy  -  Izz)  ]  ( 2-12b) 


where  R  is  the  distance  between  the  mass  centers,  is  the 
longitude  measured  from  the  x-axis,  and  is  the  latitude. 

Kinetic  Energy.  The  kinetic  energy  of  the  satellite  is 
s imply 

T  =  ( 1 / 2 ) m. v2  .  (2-13) 

1  inert 

where  the  inertial  velocity  of  the  satellite  is  given  by: 


vrel  is  the  velocity  relative  to  the  asteroid,  Q  is  the 
angular  velocity  of  the  asteroid,  and  R  is  the  distance 
between  mass  centers  as  was  previously  defined. 

Using  the  spherical  coordinates  shown,  Eq.  (2-14) 
becomes : 

''inert  =  R  •  e0  +  R  [<  X  ♦  « )  cos0]ex  ♦  ReR  (2-15) 


Thus,  the  kinetic  energy,  T,  becomes: 

T  =  (1/2)  rr^  [R20  2  +  R2(\+Jl)2cos20  +  R2] 


(2-16) 


Having  found  expressions  for  the  potential  and  kinetic 


energy,  the  Lagrangian  can  now  be  formed.  Substituting 
Eqs.  (2-12b)  and  (2-16)  into  Eq.  (2-2)  results  in  the 
following  expression  for  the  Lagrangian: 

m.  o  .  o  o  •  7  Gm. 

L  =  —  [  R>  +  R  (  X  +  S2  )  cos  0  +  R )  +  - 

2  R 

Gml  2  2 

+  — -  (  (3cos  0  cos  x  -  1)  (I  +  I  -  I  ) 

4r3  yy  zz  xx' 

+  (3cos20  sin2\  -  1)  (Ixx  +  IZ2  -  Iyy) 

+  ( 3sin20  -  1)  (Ixx  +  Iyy  -  I2Z)  ]  (2-17a) 

Since  the  mass  of  the  satellite,  m^,  is  small  and  non¬ 
zero,  it  is  convenient  to  divide  it  out  and  work  with 
equations  on  a  per  unit  mass  basis.  Thus,  the  Lagrangian  is 
rewritten  as: 


L 


+ 


+ 


+ 


1  p  ^  p  i  p  p  §  p  Gm2 

-  [  RS  +  R  (  X+Sl)  cos^0  +  R^  ]  +  - 

2  R 

G  2  2 

— -  (  (3cos  0  cos  \  -  1)  (I  +  I  -  I  ) 

4r3  y  yy  zz  xx 

( 3cos20  sin2x  -  1)  (Ixx  +  IZ2  -  Iyy) 

( 3 s i n 20  -  1)  (Ixx  +  Iyy  -  I22)  ]  ( 2-1 7b) 
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The  last  task  before  substituting  into  Eq.  (2-1)  to 
find  the  Hamiltonian  is  to  find  expressions  for  the 
generalized  velocities,  q^.  This  is  easily  accomplished  by 


using  Lagrange's  Equations 


P.  =  — T  =  K  (  X  +  n  )  COS  % 

A  o  X 


Pp  =  —  =  R  * 


P  R  =  —  =  R 

R  SR 


and  rearranging  to  find: 


X  *  - —  0 

R^cos^p 


•  p0 
0  =  - 


R  =  PR 


(2-18) 


(2-19) 


w 


111 


The  Lagrangian  can  be  rewritten  using  Eqs.  (2-19): 


1  R2p2 

<i>  2 

L  =  -  - —  +  R 


p\  \ 2  2  2  Gm2 

- - - )  COS  9  +  +  - 

R2cos20/  r  r 


2 

cos  X 

-  1)  (I 

yy 

+  Xzz  -  Xxx> 

-  1) 

(i  +  i 

'  XX  zz 

-  Iyy) 

<IXX  + 

Iyy  ~  Izz) 

] 

(2-20) 

Finally,  Eqs.  (2-19)  and  (2-20)  can  be  substituted  into 
Eq.  (2-1)  to  yield  the  Hamiltonian: 


2R2  2R2cos2$  2 


-  S)P, 


- -  [  (3cos*  cos2x  -  1)  (Ivv  +  I 

4r3  yy 

+  (3cos20  sin2x  -  1)  (Ixx  +  Izz  -  Iyy) 
+  <3sin2*  -  1)  (Ixx  +  Iyy  -  IZ2)  ] 


Xxx> 


(2-21) 


Recalling  Hamiltonian  mechanics  and  taking  the  appro¬ 
priate  partial  derivatives,  Hamilton's  canonical  equations 
of  motion  can  be  found: 


( 2-22b) 


0  = 


rr 


R  = 


(2-22 c) 


I  )  cos  o  cosx  sin\ 

yy 


( 2-22d) 


-P.  sin0  3G 

p  _  -  _  — 

©23 

R  cos  9  2R 


2  2 

7  [  (cos  X  -  sin  \ )  (I  -I  ) 
3  yy  XX 


+  (2IZZ  ■  rxx  ■  Iyy)  1  COS*  sin0 


( 2-22e ) 


•  — Gm_  3G  2  2 

P„  =  — r—  -  — t  [  (3cos  0  cos  X  -  1)  (I  +  I 
R  r2  4R'>  yy  zz 

*  (3cos2»sin2X  -  1)  (Ixx  *  Izz  -  Iyy) 

*  (3sin2,  -  1)  Uxx  *  Iyy  -  Izz)  ] 


3  2  3 

R  cos  0  R 


I 

XX 


) 


(  2  —  2  2  f ) 


Method  of  Solution 

The  problem  to  be  tackled  next  is  that  of  integrating 
the  equations  of  motion  to  find  orbits  that  are  periodic  in 
time.  This  is  done  by  choosing  the  initial  conditions  that 


will  cause  the  orbit  to  return  to  the  same  point  with  the 
same  velocity  after  the  given  period.  A  priori,  one  does 


not  know  these  initial  conditions,  so  a  method  must  be 
derived  to  find  them. 

One  logical  method  to  go  about  finding  these  proper 
initial  conditions  is  outlined  as  follows: 

1)  Select  an  orbital  period. 

2)  Guess  a  set  of  initial  conditions. 

3)  Integrate  the  equations  of  motion. 

4)  If  the  orbit  does  not  return  to  the  same  starting 

point  with  the  same  velocity,  then  vary  the 
initial  conditions  and  repeat  step  3. 

Step  1  is  simply  a  matter  of  choice.  For  Step  2, 
simple  circular  orbital  motion  is  assumed  to  give  a  fairly 
good  first  guess  at  the  initial  conditions.  The  integration 
in  Step  3  is  to  be  performed  by  Haming,  a  fourth-order 
predictor-corrector  algorithm  (6).  Thus,  the  problem  is 
reduced  to  finding  an  efficient  means  of  completing  Step  4. 

If  the  coordinates  and  momenta  are  assembled  into  the 
vector  x(t),  then  a  general  equation  of  motion  can  be 
written : 

x(t)  =  f  (x(t) ,t]  (2-23) 


Assuming  that  a  reference  solution,  xQ(t),  has  been  found,  a 
nearby  orbit  can  be  defined  as 

x  =  xQ(t)  +  <5x(t)  (2-24) 

where  6x(t)  is  a  small  displacement  in  the  reference 
orbit.  To  a  first  order  approximation,  <5x(t)  can  be  shown 
to  be  related  to  a  displacement  at  t  =  by 

x(t)  =  *(t,t]_)  6x(t1)  (2-25) 

3  x  ( t) 

Where  4(t,t.)  =  — - -  is  the  state  transition 

3  x  ( t.  ) 

xo(t) 

matrix  (6:128-129).  It  can  also  be  shown  that  the  state 
transition  matrix  can  be  found  by  solving  the  differential 
equation 

4»(  t ,  t  x  )  =  A(t)  *(  t ,  tx  )  (2-26) 

with  the  boundary  condition  that  <$(t^,t^)  =  I,  the  identity 
matrix  (6:130-131).  Aft)  is  the  variational  matrix  and  is 
defined  as: 

3  f 

A ( t)  =  —  (2-27) 

3  x 

xo(t) 

(For  reference,  the  individual  terms  of  the  matrix  are 
listed  in  Appendix  B.)  Note  that,  since  Aft)  is  evaluated 


along  the  reference  orbit,  Eq.  (2-26)  can  be  integrated 
concurrently  with  the  integration  of  Eq.  (2-23). 


Now  that  a  reference  orbit  and  the  state  transition 
matrix  have  been  found,  it  remains  to  use  this  information 
to  efficiently  vary  the  initial  conditions  to  produce  a 
closed  orbit  (Step  4).  This  is  a  relatively  simple  proce¬ 
dure  and  will  now  be  derived. 

Figure  4  shows  three  typical  orbits  starting  at  X  =  tr/2 
and  continuing  for  one-half  period.  If  these  are  truly 
periodic  orbits,  then  there  are  certain  conditions  that  must 
hold  due  to  symmetry.  These  are  clearly 

X  (t2)  =  [  X  (tj_)  +  nr] 

R29  =  “  -Vti>  (2'28) 
R  -  P„(t2)  -  PR  (tj_)  =  0 


where  X(t^)  =  ir  /  2.  These  can  be  arranged  and  placed  in  a 
vector : 


X  ( 1 2 )  -  [  X  ( t  L )  +  nr] 

Vfc2'  +  P«(tl> 
W  -  pR(ti> 


=  0  (2-29) 


If  this  vector,  G,  is  evaluated  on  a  reference  path,  x^ft), 
that  is  not  a  closed  orbit-  then  it  will  not  be  equal  to 
zero;  instead,  it  will  be  equal  to  some  error,  e.  (If  it  is 
nearly  closed,  e  will  be  small.) 


Figure  4:  Typical  Orbits 


& 


For  orbits  varied  slightly  from  the  reference  orbit, 
an  expansion  of  G  about  the  reference  orbit  can  be  written. 
To  a  first  order  approximation,  this  is  (5:6): 


3G 


G  (xQ(t2)  +  5x(t,),  t7)]  =  G  [xn(t?),  t-J  +  — 


'2''  2 


o'  2' '  2' 


3x 


<5  x  ( t2 ) 


W 


(2-30a) 


More  conveniently,  this  is  written  as 
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( 2-30b) 


G  [ x  ( 1 2 )  ,  t2J  =  G  [ XQ ( 1 2 ) ,  t2)  +  B  6x(t2) 

BG 

where  B  =  —  and  can  be  easily  be  shown  to  be: 

dx 

vv 

0  0  0  0 
0  0  0  1 
0  0  0  0 

Recall  that  on  the  reference  path: 

G  (xo(t2) ,  t2J  =  e  (2-32) 

And,  if  the  varied  orbit  is  a  closed  orbit,  then: 

G  [ x ( t 2 ) ,  t2)  =  0  (2-33) 

Thus,  after  substituting  Eqs.  (2-32)  and  (2-33)  into 
Eq.  ( 2-30b) : 

B  6  x ( 1 2 )  =  -e  (2-34) 

Recall  that  the  problem  is  to  find  the  variation 
at  t  =  tj_  that  makes  the  varied  orbit  a  closed  orbit. 
Eq.  (2-21)  relates  variations  at  t  =  t2  to  variations  at 
t  =  t^.  Thus,  Eq.  (2-34)  can  be  written: 

B  6x(t2)  =  B  ♦(t2,t1)  fixity)  =  -e 


(2-35) 


m 


The  product  B  is  a  3  x  6  matrix,  so  it  cannot 
simply  be  inverted  to  find  6x(t^).  However,  examining  this 
product  in  more  detail  reveals: 


& 


B  *(t2ft1)  6x(t1)  =  B  *(t2,t1) 


6X( tx) 
6  Q(t1) 
6R(t1) 

6P,  (tj) 
6PR(tl) 


(2-361 


Recalling  Eqs.  (2-28),  the  initial  conditions  Xft^), 
P^ft^),  and  PR(tj.)  were  set  and  not  allowed  to  vary;  there¬ 


fore,  their  corresponding  variations  at  t  =  tj_  are  zero. 


This  allows  the  first,  fifth,  and  sixth  columns  of  the 
product  to  be  eliminated  from  consideration: 


[B  $(t2,t^)  6x(t^)]224  =  [B  $  ( ^2  ,  t  ^ )  ]  2  2  4 


6  X  ( tx ) 
60  ( tx) 

5P\(ti» 


(2-37) 


Using  this  reduced  form,  Eq.  (2-35)  can  now  be  written: 


(B  ♦  ( t2  ,  ti )  1  234 


6X(tL) 


60  (tx) 


6px<t1) 


=  -e 


234 


(G[x(t2) ,t2] )234 


'.2-38) 
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The  reduced  product  is  a  square  matrix  and  can  be  inverted 
to  solve  for  6x( t^): 


ax(t1» 

«^(t1> 

6P\(V 


-l 


-  -[B  *(t2#t1)]234  EG[x(t2) ft2] )234  (2-39) 


Using  this  solution  in  Eq.  (2-24),  the  initial 
conditions  can  be  corrected.  Note  that,  since  first  order 
approximations  have  been  made  at  various  steps,  this 
correction  may  not  result  in  a  closed  orbit.  However,  if 
the  initial  reference  orbit  is  close  enough,  then  the 
correction  will  produce  a  solution  that  is  closer  to  a 
closed  orbit.  Thus,  simple  iteration  is  all  that  is  needed 
to  complete  Step  4. 

For  reference  in  later  chapters,  it  is  important  to 
note  that  the  G  vector  is  all  that  determines  the  type  of 
orbit  that  will  be  found.  The  value  of  n  will  determine  if 
it  is  a  minor  or  a  major  orbit.  (Figure  4  shows  the  result 
of  this  selection.)  The  choice  of  P^ft^)  determines  the 
inclination  of  the  orbit.  [Selecting  (t^)  =  0  results  in 
equatorial  orbits.] 


Stability  of  Orbits 

Once  a  periodic  orbit  has  been  found,  it  is  useful  to 
know  if  it  is  stable.  Since  this  is  a  Floquet  problem,  the 
condition  necessary  to  determine  stability  is  well 
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documented  (4:264-270;  6:143-148).  For  this  reason  the 
criterion  for  stability  will  be  stated  and  used  here  without 
proof . 

The  stability  of  the  orbit  is  governed  by  the  character 
of  the  Poincare'  exponents,  X^,  as  defined  by  the 
determinant 


I  *[(t1+  r),t1]  -  exp  (  X^t)  I  |  =  0  (2-40) 

where  t  is  the  period  of  the  orbit  (6:144).  Note  that  the 
terms  or  =  exptx^T)  are  simply  the  eigenvalues  of 
*  ( (t  +  r ) , t.  ]  and  that  the  Poincare'  exponents  are  given 
by: 


Xt  =  (1/t)  loge(a.)  (2-41) 

If  the  orbit  is  stable,  then  the  X^'s  must  all  be  purely 
imaginary  (4:268;  6:142). 

Notice  that,  after  integrating  Eqs.  (2-23)  and  (2-26) 
for  an  entire  orbit,  #[(t^  +  T),t^]  has  been  found. 
Thus,  it  is  a  relatively  minor  addition  to  calculate  the 
Poincare'  exponents  and,  therefore,  determine  the  stability 
of  the  orbit. 


Verification  and  Error  Detection 

Long  derivations  and  the  use  of  a  computer  invite  the 
introduction  of  errors.  With  this  in  mind,  two  tests  were 
used  help  verify  that  the  results  were  accurate. 
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Equation  Verification.  A  simple  method  can  be  derived 
to  simultaneously  verify  that  both  the  equations  of  motion, 
f,  and  the  variation  matrix,  A(t),  have  been  entered  into 
the  computer  correctly.  This  is  done  by  recalling  the 
definition: 


df 

A  { t )  =  — 

ax 


(2-42) 


xo(t) 


Each  term  of  A(t)  can  be  be  approximated  by 


a  .  . 
13 


9f  i  #  fit  ( x  j  +  A  x  j  )  ,  t  ]  -  fjjUj  -  AXj 


)  ,t] 


axj 


2  A  x  j 


(2-43) 


where  Axj  is  a  sufficiently  small  number. 

Thus,  for  any  point,  the  terms  a^  can  be  calculated 
exactly  using  the  equations  in  Appendix  B  and  approximately 
by  using  Eq.  (2-39).  If  the  values  do  not  approach  the  same 
limit  as  Axj  gets  small,  then  an  error  exists  in  either  the 
equations  of  motion  or  in  the  variation  matrix.  (It  was 
found  that  a  typical  error  resulted  in  differences  of  at 
least  25%. )  This  method  of  numerical  differentiation  was 
used  in  verifying  the  program  used  to  generate  the  results 
in  this  paper. 

Dynamics  Ver if ication  If,  during  the  numerical  inte¬ 
gration,  the  orbit  passes  too  close  to  a  singularity  in  the 
equations  of  motion,  it  is  possible  that  the  path  will  cease 
to  be  realistic.  "Too  close"  is  an  ambiguous  limitation  on 
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the  problem,  so  it  would  be  useful  if  some  way  could  be 
found  to  detect  when  it  occurs.  Fortunately,  there  is  a 
simple  method  to  detect  when  this  (or  any  other  unforeseen 
anomaly)  causes  the  integration  of  the  problem  dynamics  to 
break  down. 

Since  the  Hamiltonian  is  not  a  function  of  time  for 
this  problem,  it  is  a  constant  for  all  points  on  any  given 
orbit.  Thus,  if  the  Hamiltonian  changes  suddenly  from  one 
point  to  the  next  along  the  path,  it  is  very  probable  that 
the  limitations  on  the  dynamics  have  been  exceeded.  This 
simple  check  was  performed  automatically  at  each  integration 


s 

h 


//I 
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III.  Equilibrium  Points 


The  equilibrium  points  will  prove  to  be  important 
starting  points  in  the  computation  of  orbits.  Thus,  it  is 
advantageous  to  solve  for  them  first.  The  procedure  for 
doing  this  is  simple  and  yields  a  closed-form  solution. 


Locations  in  State  Space 

Equilibrium  points,  by  definition,  are  points  in  state 
space  from  which  the  satellite  will  not  move  if  it  is  placed 
exactly  there.  In  equation  form,  this  condition  can  be 
written : 


x  =  0 


(3-1) 


Referring  to  Eqs.  (2-22),  it  is  clear  that  a  system  of  six 
equations  in  six  unknowns  is  produced: 


2  2 
R  COS  0 


-  n  =  o 


(3-2a) 


P0 


( 3-2b) 


PR  -  0 


( 3-2c ) 


3G 

R' 


— r  (I  -  I  )  cos  0  cos\  sin\  =  0 
,3  xx  yy 


( 3  —  2d ) 
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2 

-P^  sin0  3G  2  2 

-r - 7- - -  [  (cos  \  ~  sin  X  )  (I  -  I  ) 

r2  cos30  2r3  yy  xx 


♦  (2Izz  ”  1 xx  '  Iyy)  ]  COS0  sin0  =  0  (3-2e) 


-Gw^  3G 


2  2 

~  7  (  (3cos  a  cos  X  -  1)  (I  +  I  -  I 

2  lr,4  r  yy  ZZ  XX' 


4R 


+  (3cos20sin2\  -1  )  (I  +  I  -  I  ) 

xx  zz  yy 


+  (  3  s  i  n  0  -  1)  (I  +  I  -  I  )  ] 

xx  yy  zz 


♦  —  =  0 


R^cos^^  r3 


( 3  —  2  f ) 


After  some  simple  manipulation,  the  system  can  be 
reduced  to  give  the  following  set  of  points: 


X  ,  0  ,  R ,  ,  P^  ,  Pp 


(  X  ,0  ,  R,  Px  ,  P*  ,  PR  )  =  ( 


(  0,  0,  rl,  Rj,  0,  0  ) 


0,  r2,  r2,  0,  o) 


( 3-3a ) 


(  3  —  3  b ) 


(  X  1 0  >  R )  P^  »  '  P  p  ) 


(X»0  !  R»  By  ,  P^,  Pp 


(  yr,  0 ,  R: ,  R ^ ,  0 ,  0  ) 


(3W  \ 

(  — ,  0,  r2,  r2,  0,  0 J 


(3-30 


( 3  —  3d ) 


The  terms  R,  and  R?  represent  all  of  the  physically 


realistic  roots  of  the  two  equations  that  result  from 


substituting  the  values  for  x,  $,  P^  ,  and  PR  along  with  the 
expression  for  P^  into  Eg.  (3-2f).  These  equations  are: 


3G 


-  Gm-R?  -  (I  +  I  -  21  )  =  0 

1  2  1  ^  zz  yy  xx 


( 3-4a ) 


3G 


R^  -  Gm_R?  -  (I  +  I  -  21  )  =  0 

2  21  ^  zz  xx  yy 


( 3  —  4  b ) 


Note  that  the  number  of  equilibrium  points  depends  on 
the  number  of  physically  realistic  roots  of  Eqs.  (3-4). 
Unfortunately,  it  is  impossible,  by  simple  inspection,  to 
determine  a  general  rule  regarding  the  number  of  the 
realistic  roots.  For  this  reason,  seven  actual  asteroid¬ 
like  bodies  were  examined.  The  results  were  then  analyzed 
to  see  what,  if  any,  conclusions  could  be  drawn  about  the 
typical  nature  of  the  roots  and,  therefore,  the  radii  of  the 
equilibrium  points.  (For  reference,  the  approximate  dimen¬ 
sions,  mass,  and  rotation  rate  used  for  each  of  these  bodies 
are  given  in  Table  I  (1:10;  2:138;  7:452).] 

These  calculations  did,  indeed,  reveal  information 
about  the  general  nature  of  the  roots.  Since  every  body 
examined  exhibited  the  same  trends,  it  can  be  concluded  that 
the  results  from  these  calculations  are  typical  of  realistic 
asteroids.  Thus,  the  following  generalities  can  be  made: 


aU  iTa  l*.’ ftL'ila  iUVa 1 
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Eq.  (3-4a)  has  only  one  real  root.  This  root,  Rlr 
has  a  magnitude  that  represents  a  radius  slightly 
greater  than  that  required  for  a  synchronous  orbit 
about  an  equivalent  spherical  body. 

2.  Eq.  ( 3  -  4  b )  has  three  real  roots.  Two  of  these 
roots  represent  radii  so  small  as  to  be  inside  the 
body;  therefore,  there  is  only  one  realistic  root, 
R2f  for  this  equation  also.  R2  has  a  magnitude 
slightly  less  than  that  required  for  a  synchronous 
orbit  about  an  equivalent  spherical  body. 

3.  Because  Eqs.  (3-4)  have  only  one  realistic  root 
each,  there  are  exactly  four  equilibrium  points 
given  by  Eqs.  (  3-3  )  . 

Table  I:  Typical  Asteroid  Data 


Name 

Axes  Lengths 
a  (km)  b  (km)  c  (km) 

n  (rad/s) 

Mass 

(kg) 

Hebe 

114.0 

92.0 

92.0 

2.39 

X 

io-4 

1.15 

X 

1019 

Hektor 

170.0 

63.9 

56.5 

2.53 

X 

10-4 

7.33 

X 

1018 

Juno 

137.0 

112.0 

112.0 

2.42 

X 

IO'4 

2.05 

X 

1019 

Nysa 

40.9 

27.1 

23.0 

2.73 

X 

10"4 

3.04 

X 

1017 

Pallas 

311.0 

272.0 

272.0 

2.21 

X 

10‘4 

2.75 

X 

1020 

Psyche 

142.0 

107.0 

84.9 

4.06 

X 

IO-4 

1.54 

X 

1019 

Phobos 

13.3 

11.0 

9.2 

2.28 

X 

IO-4 

1.61 

X 

1016 

1 


Stability  of  Equilibrium  Points 

An  equilibrium  point  is  stable  if  the  eigenvalues  of 
the  variation  matrix.  A,  are  all  purely  imaginary  (4:222; 
6:138-140).  The  eigenvalues  calculated  using: 

|  A  -  Ail  |  =  0  (3-5) 

where  A  is  evaluated  at  the  equilibrium  point  in  question 
and  the  Ai's  are  the  eigenvalues  at  that  point. 

It  would  be  extremely  difficult  to  find  a  general 
solution  for  the  eigenvalues  in  Eq.  (3-5).  Fortunately,  it 
is  not  necessary  to  find  a  general  solution  to  determine  the 
stability  of  the  points.  This  is  because  the  equilibrium 
points  for  the  bodies  in  Table  I  were  substituted  into 
Eq.  (3-5)  and,  once  again,  all  seven  examples  yielded  the 
same  results.  These  results,  with  reasonable  certainty,  can 
be  considered  to  be  general  for  typical  asteroids  and  are 
stated  as  follows: 

1.  The  points  at  A=  0  and  A = k  are  unstable 
equilibrium  points.  These  points  are  described 
completely  by  Eq.  (3-3a)  and  Eq.  (3-3c). 

2.  The  points  at  A  =  w/2  and  A=  (  3  ir  )  /  2  are 


stable  equilibrium  points.  These  points  are 
described  completely  by  Eq.  (3-3b)  and  Eq.  ( 3  —  3d ) . 


Example  Asteroid  Calculations 


As  an  example,  the  equilibrium  points  for  the  ficti¬ 
tious  asteroid  of  Appendix  A  were  calculated  and  checked  for 
stability.  The  results  of  these  calculations  are  summarized 
in  the  following  sections. 

Egui librium  Point  Results.  The  first  step  in  calcu¬ 
lating  these  points  was  to  find  to  roots  of  Eqs.  (3-4). 
These  roots,  as  well  as  and  R2,  are  listed  in  Table  II 
and  Table  III,  respectively.  (Note  that  the  roots  follow 

Table  II:  Roots  of  Eq.  (4-3a)  for  Example  Asteroid 


(-5.  113402188516e-01)  +  (  8 . 484056878152e-01 ) i 
(-5. 113402188516e-01)  +  (-8 . 484056878152e-01 ) i 
(  1.850707134977e-03)  +  (  2 . 466559517092e-01 )  i 
(  1.850707134977e-03)  +  ( -2 . 466559517092e-01 ) i 
(  1.018979023433e  +  00)  +  (  0 . 000000000000e  +  00 ) i 
Rx  =  1.018979023433  LU 


Table  III:  Roots  of  Eq .  (4-3b)  for  Example  Asteroid 


(-4.953656785376e-01) 

+ 

(  8 . 741840563661e-01)  i 

(-4. 953656785376e-01) 

+ 

(-8 . 741840563661e-01) i 

(-1.694701567977e-01) 

+ 

(  0.000000000000e+00) i 

(  1.703036684758e-01) 

+ 

(  0.000000000000e+00) i 

(  9.898978453971e-01) 

+ 

(  0.0J0000000000e+00) i 

R2  =  0.9898978453971  LU 

the  generalities  stated  on  page  27.)  Using  these  values  for 
R-^  and  R2,  the  equilibrium  points  were  found  with  Eqs.  (3-3) 
and  are  shown  in  Table  IV. 

Stability  Results.  The  final  conclusions  of  the  sta¬ 
bility  calculations  are  given  in  Table  IV.  While  the  deter¬ 
mination  of  stability  was  the  main  goal  of  solving  the 
eigenvalue  problem,  it  is  interesting  to  note  the 
information  that  can  be  obtained  from  the  actual  solutions 
for  the  stable  equilibrium  points. 

The  magnitudes  of  the  eigenvalues  give  frequencies  at 
which  to  start  searching  for  oscillatory  orbits  near  these 
points,  while  the  eigenvectors  give  information  about  the 


Table  IV:  Equilibrium  Points  for  Example  Asteroid 


Point 

X 

<P 

R 

P* 

PR 

Stabi 1 ity 

1 

0 

0 

R1 

PX1 

0 

0 

Unstable 

7T 

2 

2 

0 

R2 

P*2 

0 

0 

Stable 

3 

0 

R1 

PX1 

0 

0 

Unstable 

3tt 

R_ 

P,  -s 

0 

0 

Stable 

4 

IT 

0 

2 

A  2 

Rx  =  1.018979023433 

LU 

PX1  = 

1.0386121930138  LU2/TU 

R2  =  0.989897843297 

LU 

PX  2  = 

0.9801751485797  LU2/TU 
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initial  perturbation  in  direction  and  momenta  from  these 
points  (6:140-141).  [In  later  computations  of  oscillatory 
(minor)  orbits,  it  was  found  that  the  eigenvalues  and 
eigenvectors  were  of  little  use  other  than  the  determination 
of  stability.]  The  eigenvalues  and  eigenvectors  are 
listed  v n  Table  V  for  reference. 


Table  V:  Eigenvalues  and  Eigenvectors  at  Stable 
Equilibrium  Points 


Xl,2  =  (l0.4960063637745e+00) i 
(-0.9391331508767e+00)  +  (  0 . 0000000000000e  +  00 ) i 
(  0. 000000000000 Oe +00)  +  (  0 . 0000000000000e+00) i 

(  0. 000000000000 Oe +00)  +  (10 . 2894518434642e+00 ) i 
(  0. 000000000000 Oe +00)  +  (10 . 1167656750469e+00) i 
(  0.0000000000000e+00)  +  (  0 . 0000000000000e+00 ) i 
(-0.1435699563645e+00)  +  (  0 . 0000000000000e+00 ) i 


X3/4  =  (10.8673722342857e+00) i 


(-0.8799966095281e+00)  +  (  0 . 0000000000000e+00) i 
(  0.0000000000000e+00)  +  (  0 . 0000000000000e+00 ) i 
(  0.0000000000000e+00)  +  (10 . 3953216071260e+00) i 
(  0.0000000000000e+00)  +  (10 . 6043478055486e-01) i 
(  0.0000000000000e+00)  +  (  0 . 0000000000000e+00) i 
(-0.3428909856343e+00)  +  (  0 . 0000000000000e+00 ) i 


(Table  continued  on  next  page.) 


Table  V  (Continued) : 


IV.  Minor  Orbits 


Minor  orbits  are  those  orbits  that  are  merely  oscilla¬ 
tions  about  the  stable  equilibrium  points.  This  chapter 
will  present  minor  orbits  that  lie  in  the  equatorial  plane 
for  the  example  asteroid  described  in  Appendix  A.  Note  that, 
since  the  orbits  about  the  two  equilibrium  points  are  simply 
mirror  images  of  each  other,  it  is  only  necessary  to  solve 
for  the  orbits  about  one  point.  For  convenience,  the  point 
at  X  =  w/ 2  was  used. 

As  shown  in  Chapter  II,  the  G  vector  determines  the 
characteristics  of  the  orbit  that  will  be  found  by  the 
numerical  method.  If  the  integration  is  to  be  started  at 
X  =  x/2,  the  G  vector  is: 

(4-1) 

Using  this  G  vector,  minor  orbits  can  be  found  using  the 
equations  and  methods  derived  in  Chapter  II. 

Orbits  were  found  very  near  these  equilibrium  points. 
A  typical  example  of  this  family  of  orbits  is  shown  in 
Figure  5.  The  initial  conditions  are  plotted  against  the 
frequencies  of  these  orbits  in  Figures  6  and  7.  Note  that 
below  a  frequency  of  approximately  .47  rad/TU,  the  orbits 


x  ( 1 2 )  -  (ir/2) 

p»  (t2’ 
w 


became  unstable. 


Figure  5:  Typical  Minor  Orbit  Near  Equilibrium  Point 
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Figure  6:  Initial  Radius  -vs-  Frequency 
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Figure  7:  Initial  P^  -vs-  Frequency 


Another  family  of  stable  orbits  was  found  starting  at  a 
frequency  of  approximately  .44  rad/TU.  A  typical  example  of 
this  family  is  shown  in  Figure  8.  Radii  and  momenta  data 
are  given  in  Figures  9  and  10,  respectively.  Note  that  the 
orbits  of  this  family  become  unstable  near  a  frequency 
of  .395  rad/TU. 

No  other  families  of  stable  orbits  were  found.  There 
were,  however,  other  interesting  orbits  computed.  Examples 
of  a  few  of  these  are  shown  in  Figures  11-1  ilong  with 
the  initial  conditions  of  R  and  used  to  compute  them. 
(In  all  cases,  this  radius  is  the  maximum  radius  at  which 
the  orbit  crosses  perpendicular  to  the  y-axis.)  Table  VI 
summarizes  all  of  the  minor  orbits  found. 
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Figure  8:  Typical  Minor  Orbit  Near  f  =  .44  Rad/TU 


Figure  9:  Initial  Radius  -vs-  Frequency 
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Initial  Px 
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Figure  11 


,441  Rad/TU,  R  =  0.9004120106766  LU, 
0.7992199550355  LU2/TU 


Figure  14:  f  =  .290  Rad/TU,  R  =  1.354845715264  LU 


Px  =  1.117080805423  L02/TU 


Table  VI:  Summary  of  Minor  Orbits  Computed 


f 

(Rad/TU) 

R 

(LU) 

PA 

(LU2/TU) 

Stability 

0.490 

1.039041236844 

1.000867256234 

stable 

0.480 

1.062624482299 

1.012978217714 

stable 

0.470 

1.067382477069 

1.019509104377 

stable 

0.460 

1.046834141438 

1.019299715560 

unstable 

0.450 

0.986160355569 

1.005650561067 

unstable 

0.449 

0.977985586688 

1.003324500188 

unstable 

0.441 

0.900412010677 

0.977219955036 

unstable 

0.440 

0.888168589907 

0.972473210712 

unstable 

0.440 

1.238387157667 

1.056690380877 

stable 

0.435 

1.082095823329 

0.991288299428 

stable 

0.435 

1.227741787424 

1.059177135627 

stable 

0.430 

1.222281650472 

1.061495283471 

stable 

0.425 

1.219455010081 

1.063773482440 

stable 

0.420 

1.218033152644 

1.066023461238 

stable 

0.415 

1.217344653087 

1.068234383201 

stable 

0.410 

1.216985353084 

1.070392195973 

stable 

0.405 

1.216688740817 

1.072483329113 

stable 

0.400 

1.216261042734 

1.074494680583 

stable 

0.395 

1.215545133971 

1.076412660583 

stable 

0.394 

1.215354721541 

1.076783747396 

unstable 

(Table  continued  on  next  page) 


Table  VI  (Continued) 


S3 


f 

(Rad/TU) 

R 

(LU) 

PX 

(LU2/TU) 

Stability 

0.392 

1.214917630362 

1.077512354812 

unstable 

0.390 

1.214397843252 

1.078221895716 

unstable 

0.360 

1.187912621470 

1.085107385175 

unstable 

0.340 

1.117341867891 

1.077631737758 

unstable 

0.330 

1.038614678848 

1.059840481020 

unstable 

0.325 

0.984840557323 

1.044125555788 

unstable 

0.321 

1.019464633106 

0.896078589786 

unstable 

0.320 

1.029113207899 

0.897003545851 

unstable 

0.319 

1.039006180956 

0.897988374479 

unstable 

0.315 

1.081426148767 

0.902665051849 

unstable 

0.310 

1.130850261552 

0.909059307428 

unstable 

0.307 

1.195416827916 

0.918602755456 

unstable 

0.305 

1.239528500246 

0.925533794275 

unstable 

0.302 

1.350055576137 

0.942682054990 

unstable 

0.300 

1.390296759434 

1.113911631307 

stable 

0.290 

1.354845715264 

1 . 117080805423 

stable 

0.280 

1.326851733481 

1.118784303008 

unstable 

0.270 

1.296053037078 

1.119238555695 

unstable 

Major  orbits  are  those  orbits  that  completely  encircle 
the  body.  This  chapter  will  present  major  orbits  that  lie 


in  the  equatorial  plane  for  the  example  asteroid. 

Unlike  with  the  minor  orbits,  there  is  not  just  one  G 
vector  for  the  major  orbits.  The  reason  for  this  becomes 
clear  when  it  is  noted  that  the  period  of  these  orbits  is 
actually  the  synodic  period  as  defined  as 


2k 

?syn  =  |  n  -  f  j 


(5-1) 


where  S2  is  the  rotational  frequency  of  the  asteroid  and  f 
is  the  frequency  of  the  orbit.  When  f  >  n,  the  satellite 
"outruns"  the  asteroid,  implying  that  the  mean  radius  the 
orbit  is  inside  the  radius  for  synchronous  orbit. 
Similarly,  f  <  Q  implies  that  the  mean  radius  lies  outside 
the  synchronous  radius. 

Thus,  referring  back  to  Figure  3  the  G  vector  for  each 
case  can  be  easily  assembled.  If  the  integration  is  to  be 
started  at  X  =  w/2,  then  these  vectors  are: 


Xi 


X(tJ  -  (3*72) 


2) 


G  = 


X(t2)  +  (TT/2) 

W 

PR(t2} 


for  f  <  Q 


(5-2b) 


Using  these  G  vectors  and  the  method  of  solution 
derived  in  Chapter  II,  a  wide  range  of  major  orbits  were 
found.  Due  to  numerical  instability,  solutions  could  not  be 
found  in  the  range  .6  <  f  <  1.45  rad/TU.  However,  outside 
this  range,  orbits  were  computed  and  are  presented  in  the 
following  sections. 


Major  Orbits  With  f  >  fl 

The  orbits  with  a  frequency  greater  than  exhibited 
several  interesting  character istics.  The  first  is  that  they 
were  all  stable  orbits.  The  second  characteristic  is  the 
fact  that  the  paths  change  shapes  considerably  over  the 
range  of  frequencies.  Because  of  this  second  character¬ 
istic,  these  orbits  are  grouped  by  frequency  for 
presentation. 

5.2  >  f  >  1.97  Rad/TU .  An  orbit  with  a  frequency  of 
approximately  5.2  rad/TU  just  clears  the  asteroid,  so  this 
is  upper  limit  on  the  frequencies  investigated.  This  orbit 
is  shown  in  Figure  15  Figures  16  -  18  show  the  development 
of  the  family  of  orbits  as  the  frequency  decreases  (and  the 
radius  increases).  At  a  frequency  of  about  2.1  rad/TU,  the 
orbits  begin  to  rapidly  change  in  appearance. 
Figures  19  -  23  show  this  trend. 


Note  that  in  going  from  a 


frequency  of  5.2  rad/TU  to  1.97  rad/TU,  the  orbits  went 
from  being  elongated  along  the  x-axis  to  being  elongated 
along  the  y-axis. 

1.97  >  f  >  1.70  Rad/TU.  No  physically  realistic  orbits 
could  be  found  in  this  range.  All  orbits  found  had  a  shape 
similar  to  that  of  f  =  1.97  rad/TU,  with  a  radius  at  X  =  t 
so  small  as  to  cut  into  the  body. 


Figure  15:  f  =  5.2  Rad/TU,  R  =  0.3285791620856  LU, 

Px  =  0.5223527976000  LU2/TU 


Figure  16:  f  =  4.5  Rad/TU,  R  =  0.3660469358352  LU, 
Px  =  0.5586415882509  LU2/TU 


e 


Figure  17:  f  =  3.5  Rad/TU,  R  =  0.4368881991385  LUr 
P.  =  0.6178884298666  LU2/TU 


Figure  19:  f  =  2.1  Rad/TU,  R  =  0.7244296006906  LU, 
Px  =  0.7315590587046  LU2/TU 


Figure  21:  f  =  2.02  Rad/TU,  R  =  0.8294374811680  LU, 
Px  =  0.7182213754289  LU2/TU 


& 


Figure  22:  f  =  2.00  Rad/TU,  R  =  0.8676722280846  LU, 
P*  =  0.7089389646567  LU2/TU 


47 


A  .-.  A' 


Figure  23:  f  =  1.97  Rad/TU,  R  =  0.9325366690062  LUr 
Px  =  0. 6892195230246  LU2/TU 


1.70  >  f  >  1.45  Rad/TU.  In  this  range,  the  orbits  take 
on  a  simple  elliptical  shape,  as  shown  in  Figures  24  and  25. 
It  is  interesting  to  note  that  these  orbits  are,  once  again, 
elongated  along  the  x-axis. 

1.45  >  f  >  0  Rad/TU.  Due  to  these  orbits  passing  close 
to  the  stable  equilibrium  points,  the  numerical  method  used 
became  unstable.  Instead  of  converging  on  a  solution,  the 
iteration  scheme  diverged.  This  was  most  likely  due  to  a 
singularity  in  the  state  transition  matrix  caused  by  two  or 
more  very  close  solutions.  As  a  result,  no  orbits  in  this 
range  were  found. 


Figure  24:  f  =  1.70  Rad/TU,  R  =  0.5775263323319  LUf 
Px  =  0.7793806659997  LU2/TU 


$ 


Figure  25: 


f  *  1.45  Rad/TU,  R  =  0.6914873259638  LU, 
Px  =  0.8237375426793  LU2/TU 
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Major  Orbits  With  f  <  Q 

Unlike  the  inner  orbits,  these  have  almost  no  variation 
in  shape  with  frequency.  All  but  one  orbit  (f  =  .5  rad/sec) 
computed  was  stable. 

0  >  f  >  0.6  Rad/TU.  Once  again,  attempted  computation 
of  orbits  close  to  the  synchronous  radius  caused  the 
numerical  method  to  become  unstable  and  diverge.  Thus,  no 
closed  paths  in  this  range  were  found. 

0 . 6  £  Rad/TU.  As  the  frequency  decreases  (and  the 
radius  increases),  these  paths  quickly  lose  their  elongation 
and  become  almost  perfectly  circular.  This  is  shown  in 
Figures  26  and  27.  It  is  interesting  to  note  that  the  only 
unstable  major  orbit  found  occurred  at  f  =  0.5  rad/TU. 


Figure  26:  f  =  0.6  Rad/TU,  R  =  1.328869494267  LU, 
Px  =  1.198300228534  LU2/TU 
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Figure  27:  f  =  0.5  Rad/TU,  R  =  1.574940235844  LU, 
Px  =  1.271992801430  LU2/TU 


Remarks  on  Major  Orbits 

Figures  28  and  29  graphically  summarize  all  of  the 
major  orbits  computed  in  this  study.  These  computations 
reveal  several  interesting  features  of  the  major  orbits. 
These  are: 

1.  No  physically  realistic  major  with  frequencies  in 
the  range  1.97  >  f  >  1.70  rad/TU  were  found. 

2.  Interesting  phenomena  occur  near  frequencies  of 
2.0  n,  12,  and  0.512  rad/TU.  At  frequencies  very 
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asanas 


close  to  f 


2.0 fl  rad/TU,  the  orbits  begin  to 


rapidly  change  their  shape.  (See  Figures  22  -  25.) 
Near  f  =  n  rad/TU,  the  method  of  solution  became 
unstable  and  diverged.  The  only  unstable  major 
orbit  found  was  at  a  frequency  very  close  to 
f  =  0.5  SI  rad/TU. 


At  frequencies  lower  than  f  =  .4  rad/TU,  the  orbits 
are  essentially  circular.  Below  this,  the 
satellite  behaves  as  if  the  asteroid  were  a 
spherical  body. 


RADIUS  (LU) 


1.25 


1 .00 


0.75 


0.50 


0.25 


1.0  2.0  3.0  4.0  E 

FREQUENCY  ( RAD/TU) 
□FAMILY  1  AFAMILY  2  OFAMILY  3 


Figure  28:  Initial  Radius  -vs-  Frequency 
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Figure  29:  Initial  -vs-  Frequency 
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Conclusion  and  Recommendations 


The  key  results  of  this  study  have  been  stated  in  each 
chapter  already,  so  they  will  not  be  repeated  in  detail 
here.  To  summarize  the  results,  each  step  of  the  problem 
solution  will  be  looked  at  from  the  aspect  of  possible 
future  research. 

Dynamics  and  Method  of  Solution 

The  equations  of  motions  were  derived  using  a  truncated 
power  series  expansion  for  the  gravity  potential.  A  future 
study  could  retain  more  terms  in  the  series,  but  this  is  not 
recommended.  The  inaccuracy  introduced  by  ignoring  surface 
features  (craters,  etc.)  and  assuming  a  triaxial  shape  for 
the  asteroid  is  probably  greater  than  that  due  to  the  trun¬ 
cated  terms. 

If  major  orbits  near  the  synchronous  radius  are  to  be 
investigated  in  future  studies,  another  method  of  solution 
needs  to  be  found.  Instead  of  setting  a  period  and  iter¬ 
ating  on  0,  R,  and  ,  it  might  be  advantageous  to  set  R 
and  iterate  on  0,  ,  and  the  period. 

Egui 1 ibrium  Points 

The  equilibrium  points  were  fully  investigated.  No 


further  study  is  recommended. 


Minor  Orbits 


The  study  of  the  minor  orbits  produced  a  host  of 
seemingly  unrelated  paths.  These  orbits  should  be 
investigated  in  more  depth  to  see  if  any  relationships  can 
be  found.  Additionally,  inclined  orbits  should  be 
investigated. 


Major  Orbits 


If  a  new  method  cf  solution  that  remains  stable  near 
the  synchronous  radius  is  developed,  then  the  major  orbits 
in  this  region  should  be  investigated.  Once  again,  inclined 


orbits  should  also  be  investigated. 


Related  Problem  for  Study 

An  interesting  extention  of  this  research  would  be  to 
include  the  gravitational  effects  of  a  nearby  planet.  Thus, 
orbits  about  such  bodies  as  Phobos  and  Deimos  (Martian 
moons)  could  be  investigated. 


Conclusion 

This  study  found  the  governing  equations  for  the  motion 
of  a  small  satellite  about  a  rotating  asteroid.  Methods 
were  then  employed  to  solve  these  in  such  a  way  as  to  find 
the  equilibrium  points,  minor  orbits  about  the  stable  equi¬ 
librium  points,  and  major  orbits  for  a  fictitious  asteroid. 
All  orbits  found  were  in  the  equatorial  plane  of  the 


asteroid . 


Appendix  A:  Fictitious  Asteroid  Characteristics 


Any  real  asteroid  could  have  been  selected  to  use  as  an 
example,  but  Hektor  was  originally  selected  due  to  the  fact 
that  it  is  one  of  the  more  ellipsoidal  bodies.  It  was 
found,  however,  that  by  using  a  slower  rotation  rate,  the 
dynamics  produced  more  extreme  behavior.  For  this  reason, 
a  fictitious  asteroid  was  invented  with  the  dimensions  of 
Hektor  and  a  slower  rotation  rate.  The  physical  character¬ 
istics  of  this  body  are  given  in  Table  VII. 


Table  VII: 

Fictitious  Asteroid  Data 

(Standard  Units) 

a  (km) 

b  (km) 

c  (km) 

Q  (rad/s) 

Mass  (kg) 

170.0 

63.9 

56.5 

2 T  x  10'5 

7.33  x  1018 

For  numerical  reasons,  it  is  desirable  to  convert  to  a 
canonical  system  of  units  that  allows  the  integration  to  be 
performed  with  numbers  of  approximately  the  same  order  of 
magnitude.  Appropriate  length  units  (LUs),  time  units 
(TUs) ,  and  mass  units  (MUs)  can  easily  be  found  to  accom¬ 
plish  this.  If  one  length  unit  is  defined  as  the  radius  of 
a  synchronous  circular  orbit  about  an  equivalent  spherical 
body  and  if  it  is  desired  to  set  Gm2  =  1  LU^/TU^,  then  the 


following  conversion  factors  for  length,  time,  and  mass 
units  can  be  found: 


1  LU  =  467.8  km 
1  TU  =  1.592  x  104  sec 
1  MU  =  7.33  x  1018  kg 


( A— 1 ) 


After  converting  the  the  values  in  Table  VII  to  these 
units,  the  moments  of  inertia  can  also  be  calculated.  The 
physical  data  in  these  new  units  are  given  in  Table  VIII. 
These  values  were  considered  to  be  exact  and  were  used  in 
all  of  the  example  calculations. 

Table  VIII:  Fictitious  Asteroid  Data  (Canonical  Units) 


a  =  0.341  LU 
b  =  0.128  LU 
c  =  0.113  LU 


xx  ~  5.86 

X 

10“3 

lu2-mu 

yy  =  2.58 

X 

o 

1 

PO 

LU2’MU 

ZZ  =  2-65 

X 

10“2 

lu2*mu 

Q  =  .  3184  ir  rad/TU 
G  =  1  MU*LU3/TU2 


Appendix  B:  Variational  Matrix 


Terms  of  A  Matrix: 


a  =0 
11 


2  s  i  n  0 
12  R2cos30 


13  R3cos2« 


14  R2COS20 
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21 
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23  R3 


a  =0 
26 
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a36  =  1 


~ r  (i  -  i  )  (cos  \ 
r3  xx  yy 


2  2 
sm  A  )  cos  0 


6G 
7  (I 
r3  yy 


I  )  cos0  sin0  cosA  sinA 

XX 


a.,  =  —  (I  -  i  )  cos  0  cosX  sinX 
43  R4  yy  xx 


a44  =  0 


a45  =  0 


a46  =  0 


a_,  =  —  (I  -  I  )  cos0  sin0  cosX  sinX 
51  r3  yy  XX 
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52  R2  cos20  R2cos40  2R3  ^ 


+  ( 2 1 2 z  -  Ixx  -  Iyy)  ]  (cos 20  -  sin20 ) 


2P^sin0  9G  2  2 

a__  =  — ^ - x-  +  — t  [  (cos  X  -  sin  X)  (I  -  I  ) 

53  r3cos30  2r4  yy  xx 


+  (2Izz  "  Xxx  "  Iyy)  1  cos^  sin0 


2  P^  sin  0 


2  3 

R  cos  0 
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aci  =  —7  (I  -  I  )  cosX  sinx  cos  0 
61  R4  yy  XX 


2P^sin0  9G  22 

a,--  =  —7 - , - — 7  [  (I,,  -  I  )  (sin  X  -  cos  X) 
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